
## Branch and bound for H1
rm(list=ls())
library(ctgt)
library(stringr)
library(xlsx)


load("gmincmax1.Rdata")
level = tc1$level
tmin= tc1$tmin
crt=tc1$crt
cmax = tc1$cmax
levs = tc1$levs
namebank = c("1", "12", "13", "14", "15","123", "124", "125", "134", "135", "145", "1234", "1235", "1245", "1345", "12345")
h3da_1 = as.data.frame(cbind(level,tmin,crt),col.names = c("level","tmin","crt"))
h3da_2 = as.data.frame( cbind(levs,cmax),col.names = c("levs","cmax"))
line_types <- c("l1" =1,"l2" =2, "l3"=3)
point_types <- c("p1"=1, "p2"=2,"p3"=16,"p4" = 17)
library(ggplot2)
library("gridExtra")
library(gtable)
library(grid)
library(ggrepel)


h0 = ggplot(h3da_1, aes(level, tmin, shape = "p1")) + theme_classic() + xlim(level[1], level[16]+1) + ylim(0,tmin[16]+30 ) + 
  geom_point()  + 
  labs(title=expression(R%subseteq%S%subseteq%F) ) +ylab(NULL)+xlab(NULL)+ 
  theme(plot.title = element_text(hjust = 0.5, face="bold",size=16),legend.position = "none")+
  #geom_text_repel(aes(level, tmin),label=namebank,vjust = -0.5, nudge_y = 9,size= 3,segment.color = NA) +
  geom_text_repel(aes(level, tmin),label=namebank, nudge_y = 5,size= 3,segment.color = NA) +
  geom_point(aes(level, crt, shape = "p2")) + 
  #geom_text_repel(aes(level, crt),label=namebank, nudge_y = -10,size= 3,segment.color = NA)+  
  #geom_text(aes(level, crt),label=namebank,vjust = -0.5, nudge_y = -39,size= 3)+
  geom_line(aes(levs, cmax,linetype="l2"),data = h3da_2) + 
  geom_line(aes(level, tmin,linetype= "l1"), data = h3da_1[c(1,4,11,14,16),]) +
  #geom_line(aes(level, crt,linetype= "l3"), data = h3da_1[c(1,4,11,14,16),]) +
  geom_point(aes(level, tmin,shape = "p3"), data = h3da_1[c(1,4,11,14,16),]) +
  geom_point(aes(level, crt,shape = "p4"), data = h3da_1[c(1,4,11,14,16),]) +
  #labs(shape="point type", linetype = "line type") +
  scale_shape_manual(values=c(1,2,16,17), labels = c(expression(g[S]), expression(c[S]),
                                                     expression(g[tilde(R)[l]]),expression(c[tilde(R)[l]])), name="points") +
  scale_linetype_manual(values=c(1,2), labels = c(expression(g[min]), expression(c[max])), name="lines")


n= 100
m = 5
X = matrix(0, n, m,byrow = T )
for ( i in 1:n){
  set.seed(1234+i)
  X[i,] =  as.vector(arima.sim(model = list(order = c(1, 0, 0), ar = 0.2), n = m) )
}

y = rbinom(n,1,0.6)
X[which(y==1),1:3] = X[which(y==1),1:3] + 0.8


#cor(X[,1:50],X[,1:50])
library(stringr)
xs = str_replace_all(paste(rep("x",m),seq(1,m,1)),fixed(" "), "")
colnames(X) = xs

sqrW = sqrt(mean(y)*(1-mean(y)) ) # sqrt of covariance of y

WIHZ = sqrW *(sweep(X,2,colMeans(X))) ## W^{1/2}*(I-H)Z
IHZ = WIHZ/sqrW ##(I-H)Z

load("gmincmax1.RData")
namebank = c("1", "12", "13", "14", "15","123", "124", "125", "134", "135", "145", "1234", "1235", "1245", "1345", "12345")


##H1
R = "x1"

##subspace1 (F\"x3", R)
namebank2 = namebank[sapply(namebank, function(x) grepl(c("3"), x) )]
level2 =  tc1$level[sapply(namebank, function(x) grepl(c("3"), x) )]
tmin2 = tc1$tmin[sapply(namebank, function(x) grepl(c("3"), x) )]
crt2 = tc1$crt[sapply(namebank, function(x) grepl(c("3"), x) )]


R1 = c(R,"x3")
laf2 = round(eigen(tcrossprod(WIHZ[,xs,drop=F]),symmetric = T,only.values = T)$values,8) # eigen(W^{1/2}*(I-H)Z%*%Z^t%*%(I-H)%*%W^{1/2}) 
las2 = round(eigen(tcrossprod(WIHZ[,R1,drop=F]),symmetric = T,only.values = T)$values,8) 

levs2 = seq(sum(las2), sum(laf2), length.out = 200)
cmax2=c()
for (i in 1:length(levs2)){
  cmax2[i] = criticalvalue(getL(laf2,las2,levs2[i]))
}

sub2_1 = as.data.frame(cbind(level2,tmin2,crt2),col.names = c("level2","tmin2","crt2"))
sub2_2 = as.data.frame( cbind(levs2,cmax2),col.names = c("levs2","cmax2"))
line_types <- c("l1" =1,"l2" =2, "l3"=3)
point_types <- c("p1"=1, "p2"=2,"p3"=16,"p4" = 17)
library(ggplot2)

h2 = ggplot(sub2_1, aes(level2, tmin2, shape = "p1")) + theme_classic() + xlim(level2[1], level2[16]+5) + ylim(0,tmin2[16]+35 ) + 
  geom_point()  + 
  labs(title=expression(("R"~union("{3}"))%subseteq%S%subseteq%F) ) +ylab(NULL)+ xlab(NULL)+ 
  theme(plot.title = element_text(hjust = 0.5, face="bold",size=16),legend.position = "none")+
  #geom_text_repel(aes(level2, tmin2),label=namebank2,vjust = -0.5, nudge_y = 9,size= 3,segment.color = NA) +
  geom_text_repel(aes(level2, tmin2),label=namebank2, nudge_y = 5,size= 3,segment.color = NA) +
  geom_point(aes(level2, crt2, shape = "p2")) + 
  #geom_text_repel(aes(level2, crt2),label=namebank2, nudge_y = -10,size= 3,segment.color = NA)+  
  #geom_text(aes(level2, crt2),label=namebank2,vjust = -0.5, nudge_y = -39,size= 3)+
  geom_line(aes(levs2, cmax2,linetype="l2"),data = sub2_2) + 
  geom_line(aes(level2, tmin2,linetype= "l1"), data = sub2_1[c(1,3,7,8),]) +
  #geom_line(aes(level, crt,linetype= "l3"), data = sub1_1[c(1,4,11,14,16),]) +
  geom_point(aes(level2, tmin2,shape = "p3"), data = sub2_1[c(1,3,7,8),]) +
  geom_point(aes(level2, crt2,shape = "p4"), data = sub2_1[c(1,3,7,8),]) +
  #labs(shape="point type", linetype = "line type") +
  scale_shape_manual(values=c(1,2,16,17), labels = c(expression(g[S]), expression(c[S]),
                                                     expression(g[tilde(R)["\u2113"]]),expression(c[tilde(R)["\u2113"]])), name="points") +
  scale_linetype_manual(values=c(1,2), labels = c(expression(g[min]), expression(c[max])), name="lines")



#subspace2 (F, c(R,"x3))
namebank1 = namebank[!sapply(namebank, function(x) grepl(c("3"), x) )]
level1 = tc1$level[!sapply(namebank, function(x) grepl(c("3"), x) )]
tmin1 = tc1$tmin[!sapply(namebank, function(x) grepl(c("3"), x) )]
crt1 = tc1$crt[!sapply(namebank, function(x) grepl(c("3"), x) )]



xs1 = setdiff(xs,"x3")
laf1 = round(eigen(tcrossprod(WIHZ[,xs1,drop=F]),symmetric = T,only.values = T)$values,8) # eigen(W^{1/2}*(I-H)Z%*%Z^t%*%(I-H)%*%W^{1/2}) 
las1 = round(eigen(tcrossprod(WIHZ[,R,drop=F]),symmetric = T,only.values = T)$values,8) 

levs1 = seq(sum(las1), sum(laf1), length.out = 200)
cmax1=c()
for (i in 1:length(levs1)){
  cmax1[i] = criticalvalue(getL(laf1,las1,levs1[i]))
}

sub1_1 = as.data.frame(cbind(level1,tmin1,crt1),col.names = c("level1","tmin1","crt1"))
sub1_2 = as.data.frame( cbind(levs1,cmax1),col.names = c("levs1","cmax1"))
line_types <- c("l1" =1,"l2" =2, "l3"=3)
point_types <- c("p1"=1, "p2"=2,"p3"=16,"p4" = 17)
library(ggplot2)



h1 = ggplot(sub1_1, aes(level1, tmin1, shape = "p1")) + theme_classic() + xlim(level1[1], level1[16]+5) + ylim(0,tmin1[16]+35 ) + 
  geom_point()  + 
  labs(title=expression(R%subseteq%S%subseteq%F~"\\ {3}")) + ylab(NULL)+ xlab(NULL)+ 
  theme(plot.title = element_text(hjust = 0.5, face="bold",size=16),legend.position = "none")+
  #geom_text_repel(aes(level1, tmin1),label=namebank1,vjust = -0.5, nudge_y = 9,size= 3,segment.color = NA) +
  geom_text_repel(aes(level1, tmin1),label=namebank1, nudge_y = 5,size= 3,segment.color = NA) +
  geom_point(aes(level1, crt1, shape = "p2")) +  
  #geom_text_repel(aes(level1, crt1),label=namebank1, nudge_y = -10,size= 3,segment.color = NA)+ 
  #geom_text(aes(level1, crt1),label=namebank1,vjust = -0.5, nudge_y = -39,size= 3)+
  geom_line(aes(levs1, cmax1,linetype="l2"),data = sub1_2) + 
  geom_line(aes(level1, tmin1,linetype= "l1"), data = sub1_1[c(1,3,7,8),]) +
  #geom_line(aes(level, crt,linetype= "l3"), data = sub1_1[c(1,4,11,14,16),]) +
  geom_point(aes(level1, tmin1,shape = "p3"), data = sub1_1[c(1,3,7,8),]) +
  geom_point(aes(level1, crt1,shape = "p4"), data = sub1_1[c(1,3,7,8),]) +
  #labs(shape="point type", linetype = "line type") +
  scale_shape_manual(values=c(1,2,16,17), labels = c(expression(g[S]), expression(c[S]),
                                                     expression(g[tilde(R)["\u2113"]]),expression(c[tilde(R)["\u2113"]])), name="points") +
  scale_linetype_manual(values=c(1,2), labels = c(expression(g[min]), expression(c[max])), name="lines")



library("gridExtra")
library(gtable)
library(grid)
# png("bb.png",width = 22, height = 18, units = 'cm', res = 1200)
# grid.arrange(arrangeGrob(ggplot()+ theme(panel.background = element_blank()), h0,ggplot()+ theme(panel.background = element_blank()),ncol=3,widths = c(1,2,1)),                          
#              arrangeGrob(ggplot()+ theme(panel.background = element_blank()),
#                          ggplot() + geom_segment(data=data.frame(x=1, y=1, vx=-0.2, vy=-0.2), mapping=aes(x=x, y=y, xend=x+vx, yend=y+vy), arrow=arrow()) +theme_void() ,
#                          ggplot()+ theme(panel.background = element_blank()),
#                          ggplot() + geom_segment(data=data.frame(x=1, y=1, vx=0.2, vy=-0.2), mapping=aes(x=x, y=y, xend=x+vx, yend=y+vy), arrow=arrow()) +theme_void(),
#                          ggplot()+ theme(panel.background = element_blank()),          
#                          ncol = 5,widths = c(2,1,2,1,2)),
#              arrangeGrob(h1, ggplot()+ theme(panel.background = element_blank()),h2, ncol = 3,widths = c(5,1,5)), # Second row with 2 plots in 2 different columns
#              nrow = 3,
#              heights=c(5,0.8,5),
#              left = textGrob("test statistics & critical values", rot = 90, vjust = 1),
#              bottom = textGrob("levels"))
# 
# dev.off()
# 
# 
# 
# png("bb_h0.png",width = 16, height = 12, units = 'cm', res = 800)
# h0
# dev.off()
# 
# png("bb_h1.png",width = 16, height = 12, units = 'cm', res = 800)
# h1
# dev.off()
# 
# png("bb_h2.png",width = 16, height = 12, units = 'cm', res = 800)
# h2
# dev.off()
# 
tiff("bb.tiff",width = 16, height = 14, units = 'cm', res = 600)
#jpeg("bb.jpg",width = 16, height = 14, units = 'cm', res = 600)
grid.arrange(h0,
             arrangeGrob(h1, ggplot()+ theme(panel.background = element_blank()),h2, ncol = 3,widths = c(5,1,5)), # Second row with 2 plots in 2 different columns
             nrow = 2,
             left = textGrob("test statistic & critical value", rot = 90, vjust = 1),
             bottom = textGrob("level"))

dev.off()

